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NOTE 

A  Three-Point  Sixth-Order  Nonuniform 
Combined  Compact  Difference  Scheme 


A  three-point  nonuniform  combined  compact  difference  (NCCD)  scheme  with 
sixth-order  accuracy  is  proposed  for  numerical  models.  The  NCCD  scheme  is  a  gen¬ 
eralization  of  the  previously  proposed  combined  compact  difference  (CCD)  scheme 
with  a  global  Hermitan  polynomial  spline  and  has  major  improved  features  such  as 
error  and  computational  (CPU)  time  reduction.  For  nonperiodic  boundaries,  addi¬ 
tional  sixth-  or  fifth-order  nonuniform  boundary  conditions  are  proposed.  The  NCCD 
scheme  with  either  sixth-  or  fifth-order  additional  boundary  conditions  can  increase 
the  accuracy  and  decrease  the  CPU  time  about  1-2  orders  of  magnitude,  compared 
to  the  CCD  scheme.  ©  1999  Academic  Press 


1.  INTRODUCTION 

Following  the  trend  toward  highly  accurate  numerical  schemes  of  partial  differential 
equations  (PDE)  led  by  many  authors  (e.g.,  Adam  [1];  Chu  and  Fan  [2];  Hirsh  [5];  Rubin 
and  Khosla  [8];  Navon  and  Riphagen  [6]),  Chu  and  Fan  [2]  proposed  a  three -point  sixth- 
order  uniform  combined  compact  difference  (CCD)  scheme  to  increase  accuracy.  This 
scheme  follows  earlier  work  on  the  use  of  second  derivatives  in  compact  differencing  (such 
as  Rubin  and  Khosla  [8]), 


1 

Y  (akfi+k  +  bkf!+k  +  ckf"+k)  =  0,  (1) 

k=-\ 

which  is  referred  as  the  Hermite  formula.  Here,  /  is  a  dependent  variable. 

The  ocean  variability  is  not  uniform.  For  example,  the  western  boundary  currents  (e.g., 
the  Gulf  Stream,  Kuroshio)  have  much  larger  variability  than  the  ocean  interior.  An  ideal 
treatment  is  to  use  a  nonuniform  scheme:  high  resolution  grids  for  high  variability  areas 
and  low  resolution  grids  for  low  variability  areas.  Goedheer  and  Potters  [4]  proposed  a 
nonuniform  grid  for  a  fourth -order  compact  difference  scheme.  Following  their  path,  we 
propose  in  this  paper  a  three -point  sixth-order  nonuniform  combined  compact  (NCCD) 
scheme  with  sixth-order  or  fifth-order  accuracy  at  both  interior  and  exterior  boundaries. 
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2.  NCCD  SCHEME 

Let  fix)  be  defined  on  the  interval,  0  <  x  <  L.  Use  a  nonuniform  grid,  0  =  xo  <  x\  <  X2  < 
■  ■  ■  <  xN  =  L  with  a  nonuniform  space  hi  =  Ax,  =  x,  —  x,_i  (i  =  1,  2, . . . ,  N).  Let  the  de¬ 
pendent  variable  fix)  at  any  grid  point  x,  and  two  neighboring  points  x,_i  and  x,+l  be  given 
by  fi,  fi-i,  and  f+\  and  let  its  derivatives  at  the  two  neighboring  points  x,_i  and  x,+i  be 
given  by  f[_x,  f!'_v  and  f'+l,  f"+1.  Let  H,(x)  be  a  local  Hermitian  polynomial  defined  on 
a  closed  interval  [x/_i,x,-+i]  with  a  spacing  of  (hj  +/z,+i),  and  determined  by  /  evaluated 
at  x,-_i,  Xj,  Xi+ 1,  and  its  derivatives  /',  and  f"  evaluated  at  x,_ i  and  x,-+i  (Fig.  la), 

H,(x)  =  f  +  (/;■_!  -  mm  +  ifi+ 1  -  +  fUhiQiti) 

+  f!+lhi<S>  4(f)  +  fUh^siH)  +  f"+lhj^),  (2) 

where  f  =  (x  —  x,  )//z,-  and 

4>y«)  =  a^  +  ^2  +  ^3  +  ^4  +  ey?5  +  ^6.  7  =  1,2 . 6,  (3) 

are  six  elements  of  the  local  Hermitian  polynomial  W,  (x)  satisfying  the  following  homo¬ 
geneous  conditions 

®;(-l)=0,  <!>;(£,■)  =  0,  <*>'■(—!)  =  0,  %(k,)  =  0, 

0"(-l)  =  0,  <t>"j{ki)  =  0,  j  =  1,2, ...  ,6, 

except  for 

4»t(-l)=l,  4>2(*/)  =  l,  <*>3(— 1)  =  1,  <*>;(*/)  =  1,  4»?(-l)=l,  <(*,-)  =  1, 


(a) 


i-1  i  i+1 


0  12  3 

FIG.  1.  Nonuniform  grid  systems:  (a)  three-point  grid  in  the  interior  domain;  (b)  three-point  grid  at  the  left 
boundary;  and  (c)  four-point  grid  at  the  left  boundary. 
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where  k,  =h,+\/h,. 
If  we  define 


l  —  l 


Sf 


=  fi-vr 


=  f! 


+1> 


i+1 


S2f 


±  =  */(*),  hdr  = + ++ 


8x2 

s2f 

8x 2 


=  fi—  I  ’ 


i-1 


+/ 

8x2 


=  f!U. 


i+ 1 


the  NCCD  scheme  is  given  by 


kf(2kj  +  5)  f8f 
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(ki  +  l)4  \SxJi+ 1 


52/ 


2(*/  +  l)3  V  5a-2  y 2(kj  +  l)3  \8x2Ji+l 
3  ($kf  +  4 kt  +  1)  f,+l  -  ft  ,  3*?(U?  +  4*,  +  5)  f  -  /J.j 


(*,-  +  l)5 


kjhj 


&  +  l)5 


(4) 


and 


6*,?  (5 -*,•-*,?)  (8f 
(kt  + 1)4  V«jc 


6(5*,? /«/ 


I— 1 


k2(3-2ki)hf82l 


ki(ki  +  l)4  V^/, 
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fi+ 1  ~  fi 
k;h : 


hi 


(5) 


Expanding  the  variable  /,  /',  and  f"  into  the  Taylor  series  at  points  i  —  1  and  i  +  1,  the 
truncation  errors  are  estimated  as  —kff^’hf/ll  for  (4)  and  —2k2 {9k2  —  11  kj  +  9)ffS)h]/8\ 
for  (5). 


3.  A  GLOBAL  HERMITAN  POLYNOMIAL  SPLINE  ON  A  NONUNIFORM  GRID 

Let  H, (x )  be  a  local  Hermitian  polynomial  defined  on  a  closed  interval  [x,_i ,++i]  with 
a  spacing  of  (hi  +  hl+\ ).  We  may  define  a  global  Hermitan  polynomial  spine  by 

Hg(x)  —  Hi(x),  0  =  X\  <  x  <  xi, 

Hg(x)  —  - Hi(x)  H - - — -Hi+ i(x),  x,  <  x  <  xl+\  (i  =  2,  3,  . . . ,  N  -  2)  (6) 

hi 

Hg(x )  =  Hn-i(x),  xn- i  <  x  <  jcjv, 


which  has  up  to  third-order  continuous  derivatives  at  [vt,  xn]. 
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4.  BOUNDARY  ALGORITHMS 

For  periodic  boundaries,  the  NCCD  scheme  automatically  provides  the  sixth-order  ac¬ 
curacy  at  the  boundaries.  But  for  nonperiodic  boundaries,  we  propose  a  sixth-order  and 
a  fifth-order  boundary  algorithm  for  solving  finite  difference  equations  (FDE)  and  a  fifth- 
order  boundary  algorithm  for  difference  calculation.  For  simplicity,  we  discuss  left  boundary 
(i  —  0)  as  an  example.  The  treatment  of  the  right  boundary  (i  —  N)  is  the  as  same  as  for  the 
left  boundary. 


4.1.  Finite  Difference  Equation  (FDE) 

We  have  a  choice  of  using  sixth-order  or  fifth-order  additional  boundary  conditions.  The 
grid  structure  is  illustrated  in  Fig.  lb. 

4.1.1.  Sixth-Order  Accuracy  Formulation 

The  sixth-order  accuracy  at  the  boundary  is  based  on  using  a  global  Hermitan  polyno¬ 
mial  spline  and  an  integration  equation  for  the  boundary  cell.  For  a  second-order  partial 
differential  equation  with  constant  coefficients, 

af'ix)  +  bf'(x)  +  c/(x)  =  s(x),  (7) 


we  propose  an  extra  boundary  condition  with  sixth-order  accuracy  as 


(s£)  0  + ' 02  (sx  )  ,  +  bl  (si)  o  +  bl  (si)  , +  Ci  /o  +  C2/l  +  C3/2  -  d •  (8) 

4.1.2.  Fifth-Order  Accuracy  Formulation 

If  we  can  bear  a  little  less  accuracy  (fifth-order)  at  the  boundary,  the  formulation  will  be 
much  simpler. 
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=  0. 


(9) 


4.2.  Finite  Difference  Calculation 

The  NCCD  can  also  be  used  to  calculate  the  high-order  finite  difference.  There  are 
2  *  (N  +  1)  unknown  variables  [(S//Sx);,  (S2//Sx2),,  i  =  0,  1,  2, . . . ,  N],  but  there  only 
are  2  *  (N  —  1)  equations  (4)  and  (5)  at  internal  nodes.  Therefore,  at  each  boundary  point 
two  additional  conditions  are  needed  to  close  the  system.  Here,  we  only  show  the  left 
boundary.  The  first  additional  boundary  condition  is  (9).  Expanding  the  variable  /  into 
the  Taylor  series  in  another  form  (Fig.  lc),  we  obtain  the  second  additional  boundary 
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condition. 
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(10) 


5.  STOMMEL  OCEAN  MODEL 

We  use  the  Stommel  ocean  model  [9]  to  compare  the  accuracy  and  CPU  between  NCCD 
and  CCD  schemes. 


5.1.  Model  Description 

Stommel  [9]  designed  an  ocean  model  to  explain  the  westward  intensification  of  wind- 
driven  ocean  currents .  Consider  a  rectangular  ocean  with  the  origin  of  a  Cartesian  coordinate 
system  at  the  southwest  corner.  The  x  and  y  axes  point  eastward  and  northward,  respectively. 
The  boundaries  of  the  ocean  are  at  x  =  0,  Lx  and  y  =  0,  Ly.  The  ocean  is  considered  as  a 
homogeneous  and  incompressible  layer  of  constant  depth  H  when  at  rest.  Stommel  derived 
an  equation  for  the  streamfunction  i // , 

/  32  32  \  (n  \ 

X^-  +  wT+l>^  =  -vim\v)  (11) 

with  the  boundary  conditions: 

*(0,  y )  =  V{LX,  y)  =  *(x,  0)  =  vp(x,  Ly)  =  0.  (12) 

Here,  the  right-hand  side  of  ( 1 1)  indicates  the  wind  forcing,  and  fi ,  y  represent  the  latitudinal 
change  of  the  Coriolis  parameter  and  the  strength  of  the  surface  wind  stress,  respectively. 
The  analytical  solution  of  (11)  with  the  boundary  conditions  (12)  is  given  by  (Fig.  2) 

^(a\x,y)  =  -y(^-^j  sm(j*-y^(peAx  +qeBx  -l).  (13) 


FIG.  2.  Streamfunction  (m2/s)  obtained  from  Stommel  ocean  model. 
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The  physical  parameters  selected  by  Stommel  [9]  are 

Lx  =  107  m,  Ly  =  2jt  x  106m,  ft  =  ^  x  10_5m_1,  y  =  ^  x  10“10m-2.  (14) 

5.2.  Two-Dimensional  Nonuniform  Grid 

We  discretize  the  domain  into  0  =  xo  <  x\  <  X2  <  ■  ■  ■  <  xM  =  Lx ;  0  =  y0  <  y  i  c  \n  <  ■  ■  ■ 
<yN  —  Ly  with  nonuniform  spacing  (Ax,-,  Ay,),  that  is, 

Xi  =  x,_i  +  Ax,-,  yj  =  y,_i  +  Ay,-,  i  =  1,2,  j  =  1,2, ...,  N,  (15) 

where  M  and  N  are  numbers  of  grid  cells  in  the  x  and  y  directions,  respectively.  A  two- 
dimensional  nonuniform  grid  (Fig.  3)  was  used  with  left-to-right  decreasing  resolution  in 
the  x-direction, 


Ax,-  =  Axoo 


i  =  1,2, 


(16) 


and  with  boundary-to-interior  decreasing  resolution  and  symmetric  at  the  middle  in  the 
y-direction, 


Ay,-  =  Ayoo 


r  (  2A] 

1  ...  . 

~  N~ 

1  —  exp  ay—  J 

VI 
— ■, 

VI 

4— < 

_y_ 

Ay,-  =  Ayw_,-+i, 


<j<N , 


(17) 


X  (1 06  m)  grid  size 


FIG.  3.  Two-dimensional  nonuniform  grid. 
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where  ax  and  are  the  e -folding  homogeneity  parameters  and  Axx  and  A are  the 
limiting  grid  size  in  the  x  and  y  directions,  respectively.  Increase  of  ax  (or  ay )  means 
the  increase  of  the  homogeneity  of  the  x-grid  (or  y-grid).  The  values  of  Axqo  and  Ay^  are 
determined  by  satisfying  the  conditions  xM  =  Lx  and  yN  —  Ly.  In  this  study,  we  used  the 
same  number  of  grid  cells  for  both  x  and  y  directions  (M  —  N ). 

We  use  the  same  alternating  direction  implicit  (ADI)  method  as  described  in  Chu  and 
Fan  [3]  to  solve  the  two-dimensional  equation  (11)  numerically.  Such  an  iterative  process 
stops  when  the  correction  at  the  iteration  k  +  1, 


r(*+D  _ 


EA/-1  Y^w-i 

i=i  l^j= 1  ^i,j 


(A  Xi  +  Ax/+i)(A  yj  +  Ayj+1) 


(Ax>  +  Ax!+i)(Ay7-  +  Ayy-+1) 


is  smaller  than  10  8.  The  numerical  solution  for  the  grid  point  (x, ,  y  f )  is  T,  ( .  The 


TABLE  I 

Comparison  of  Relative  Average  Error  (%) 


n=15  (CCD  AARE=  0.00236) 


n=20  (CCD  AARE=  0.000238) 
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discretization  error  was  represented  by  an  area-averaged  relative  error  (AARE), 


AARE  = 


EEEE1 


(A  Xi  +  Axi+l)(Ayj  +  Ayj+1) 


E«  E^TC  KA*i  +  V'V-,HAv;  +  Ayj+l) 


We  solved  (11)  numerically  with  both  NCCD  and  CCD  schemes  under  various  horizontal 
resolutions  and  recorded  the  CPU  time  (a  SUN  Sparc -20  was  used)  for  each  run.  To  test 
the  sensitivity  of  NCCD  scheme  on  the  grid  size  and  e-folding  scale,  we  computed  the  dis¬ 
cretization  error  with  different  cell  numbers  M  —  N  —  15, 20, 25, 30, 35, 40,45, 50,  different 
a; -directional  e-folding  scales  ax  =  0.1,  0.2,  0.3,  0.5,  0.7,  1,  2  and  different  y-directional  e- 
folding  scales  ay  =  0.2,  0.5,  1 .  Using  the  uniform  CCD  scheme  as  the  reference,  we  define 
the  error  and  CPU  ratios  (%)  by 

AARE(NCCD)(A,  ax,  ay)  CPUfNCCD)((V,  ax,  ay) 

e(N,  ax,  av)  =  - - — ,  X(ax,  av)  —  - - — .  (20) 

-  AARE<CCD)(A^)  •  CPU(CCD)(A) 

Thus,  the  uniform  CCD  scheme  has  a  value  of  100%  for  both  e  and  X. 


TABLE  II 

Comparison  of  Relative  CPU  Time  (%) 


n=45  (CCD  CPU  time=  1 69(s) )  n=50  (CCD  CPU  time=  274(s) ) 
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a  =0.2 
y 


a  =0.5 
y 


Y-uniform 


FIG.  4.  The  AARE-CPU  diagrams  for  NCCD  and  CCD  scheme  comparison.  Here,  the  solid  curve  indicates 
the  CCD  scheme.  Eight  symbols  represent  different  a,-values:  0.1  (o);  0.2  (A);  0.3  (□);  0.5  (o);  0.7  (+);  1.0  (x); 
2.0,  (*);  and  oo  (v)- 


Table  I  shows  the  distribution  of  s(N,  ax,  ay).  The  symbols  “X-D”  and  “Y-D”  mean  the 
x-dependency  and  v-depciidcncy,  respectively.  The  last  column  (row)  of  the  table  indicates 
the  use  of  the  uniform  grid  in  x -direction  (y -direction).  The  AARE  of  the  NCCD  scheme 
is  1  to  14%  of  the  AARE  of  the  uniform  CCD  scheme.  This  error  reduction  enhances  as 
N  increases  and  as  ax  or  ay  decreases.  We  further  notice  that  the  AARE  are  very  sensitive 
to  the  e-folding  parameter  ax  and  very  insensitive  to  ay.  Taking  A'  =  1 5  as  an  example, 
s  decreases  from  14.2%  to  8.35%  as  ax  decreases  from  2  to  0.1.  However,  s  keeps  the 
same  value  as  ay  varies.  Such  an  x—y  asymmetry  is  called  by  the  inhomogeneous  ocean 
variability.  The  analytical  solution  of  the  Stommel  ocean  model  shows  a  strong  variability 
of  T  in  the  western  boundary  only  (Fig.  2).  The  nonuniform  scheme  reduces  the  truncation 
error  drastically  when  the  fine  resolution  grid  is  applied  there. 

Table  II  shows  the  distribution  of  X(N,  ax,  ay).  The  CPU  of  the  NCCD  scheme  is  50  to 
89%  of  the  CPU  of  the  uniform  CCD  scheme.  This  CPU  saving  enhances  as  N  decreases 
and  as  ax  increases.  Taking  N  =  15,  ay  =  0.5  as  an  example,  /,  decreases  from  54.6%  to 
50.8%  as  ax  increases  from  0.1  to  2. 

We  use  the  AARE-CPU  diagram  to  verify  the  performance  of  the  NCCD  scheme.  When  ax 
and  ay  are  given,  both  AARE  and  CPU  depend  on  N.  Thus,  variation  of  N  creates  a  curve  in 
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(a) 


CPU  time  (s) 

FIG.  5.  Error  comparison  between  specific  NCCD  scheme  (with  a ,  =  0.3  and  ay  =  0.5)  and  CCD  scheme: 
(a)  CPU  time;  (b)  AARE;  and  (c)  AARE-CPU  diagram.  Here,  the  solid  curves  denote  the  CCD  scheme,  and  the 
dashed  curves  denote  the  NCCD  scheme. 


the  AARE-CPU  diagram.  We  only  use  eight  different  values  for  A  (15, 20, 25, 30, 35, 40, 45, 
50)  in  this  study.  Therefore,  we  have  eight  points  given  for  ax  and  ay  on  the  diagram.  Figure  4 
shows  the  AARE-CPU  diagram  for  four  different  values  of  <xy  (0.2, 0.5, 1 ,  and  oo  (uniform)), 
respectively.  The  solid  curve  indicates  the  CCD  scheme.  Eight  different  symbols  represent 
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the  eight  different  values  used  for  ax .  The  symbol  “V”  indicates  ax=oo  (uniform).  The  other 
seven  symbols  represent  the  seven  values  we  used  for  ax  (0.1, 0.2, 0.3, 0.5, 0.7,  1, 2).  For  the 
CCD  scheme,  ax  —  oo  and  av  =  oo.  Again,  we  see  the  insensitivity  of  the  NCCD  scheme  to 
ay.  For  a  given  CPU  time,  the  NCCD  scheme  increases  the  accuracy  by  1-2  orders  of  mag¬ 
nitudes. 

Figure  5  shows  the  comparison  between  a  specific  NCCD  scheme  (with  ax  =  0.3  and 
ay  =  0.5)  and  the  CCD  scheme.  The  solid  curves  denote  the  CCD  scheme,  and  the  dashed 
curves  denote  the  NCCD  scheme.  This  specific  NCCD  scheme  reduces  both  CPU  time 
by  20-30%  (Figure  5a)  and  AARE  around  100  times  (Fig.  5b)  versus  the  CCD  scheme. 
Such  a  reduction  enhances  as  N  increases.  Such  an  error  reduction  is  clearly  seen  in  the 
AARE-CPU  diagram  for  this  specific  NCCD  scheme  and  the  CCD  scheme  (Fig.  5c). 


6.  CONCLUSIONS 

(1)  This  study  shows  that  the  NCCD  scheme  is  a  promising  highly  accurate  method 
for  both  derivative  computation  and  FDE  solutions.  The  NCCD  scheme  has  all  the  good 
features  as  the  CCD  scheme,  such  as  three-point  and  sixth-order  accuracy.  Besides,  the 
NCCD  scheme  represents  the  high  variability  area  more  accurate  than  the  CCD  scheme  by 
using  fine  grid  there.  Our  analysis  shows  that  for  the  same  CPU  time,  the  NCCD  scheme 
has  an  error  reduction  by  1-2  orders  of  magnitude,  compared  to  the  CCD  scheme. 

(2)  To  keep  a  three-point  scheme,  an  additional  boundary  algorithm  is  needed.  We  pro¬ 
pose  both  fifth-  and  sixth-order  boundary  algorithms.  Since  the  fifth-order  formulation  is 
more  simple  than  the  six-order  formulation  and  easy  to  use,  especially  in  two-dimensional 
problems. 

(3)  Use  of  the  NCCD  scheme  makes  the  convergence  fast.  In  other  words,  it  needs  less 
iterations  and  saves  CPU  time. 

(4)  We  found  a  global  Hermitan  polynomial  spine  in  this  study.  Thus,  the  NCCD  scheme 
is  easily  used  in  integral  methods  as  well  as  in  stagged  grids. 
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